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Abstract 

Recently we developed a diagonal homotopy method to compute a 
numerical representation of all positive dimensional components in the 
intersection of two irreducible algebraic sets. In this paper, we rewrite 
this diagonal homotopy in intrinsic coordinates, which reduces the 
number of variables, typically in half. This has the potential to save a 
significant amount of computation, especially in the iterative solving 
portion of the homotopy path tracker. Three numerical experiments 
all show a speedup of about a factor two. 
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Our goal is to compute the irreducible decomposition of A fl B C C k , 
where A and B are irreducible algebraic sets. In particular, suppose that 



• A is an irreducible component of the solution set of a polynomial system 

= defined on C fc , and similarly 

• B is an irreducible component of the solution set of a polynomial system 
Jb(u) = defined on C k . 

This includes the important special case when /a and fs are the same system, 
but A and B are distinct irreducible components. 

Casting this problem into the framework of numerical algebraic geom- 
etry, we assume that all components are represented as witness sets. For 
an irreducible component A C C fc of dimension dim(A) and degree deg(A), 
a witness set consists of a generic k — dim (A) dimensional linear subspace 
L C C k and the deg(A) points of intersection An L. We assume that at 
the outset we are given such sets for A and B, and our goal is to compute 
witness sets for the irreducible components of An B. The intersection may 
break into several such components, and the components may have various 
dimensions. Our methods proceed in two phases: we first find a witness su- 
perset guaranteed to contain witness points for all the components, then we 
break this set into its irreducible components. We recently reported on an 
algorithm ^3], herein called the extrinsic 1 homotopy method, for computing 
a witness superset for An B. This can then be decomposed into irreducible 
components using the methods in ^3] and its references. 

Abstracting away the details, which are discussed more fully in Jfl the 
extrinsic method consists of a cascade of homotopies in unknowns x G C N 
and path parameter t 6 [0,1], each of the form 



H(x,t) 



t(Px + p) + (l-t)(Qx + q) 



(1) 



where / : — > C m is a system of polynomial equations, P, Q are (N — 
m) x iV full-rank matrices, and are column vectors. There is 

a homotopy of this form for each dimension where An B could have one or 
more solution components. We know solution values for x at t — 1 and wish 
to track solution paths x(t) implicitly defined by ((TJ) as t — » to get x(0). 

1 The terminology extrinsic/intrinsic is in analogy with the homotopies of 0]. 
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At any specific value of t, this looks like 



H(x,t) 



R(t)x + r(t) 







(2) 



where R = tP + (1 —t)Q and r = tp+ (1 —t)q. The homotopy is constructed 
such that we are assured that R(t) is full rank for all t G [0, 1]. Thus, the 
linear subspace of solutions of R(t)x + r(t) = can be parameterized by 
u G C m in the form 



where x p (t) is any particular solution and R (t) is the right null space of 
R(t), that is, R 1 - is a full-rank N x m matrix with RR 1 - = 0. We may 
restrict H to this linear subspace to obtain 



where we have dropped the linear equations because by construction, they 
are identically zero for all t. We refer to this as the intrinsic form of the 
equations. 

The problem with (@J) is that it requires computing R x and x v at each 
new value of t as we follow the homotopy paths. Because of this, H(x) offers 
little, if any, computational advantage over the extrinsic H(x). 

Although not generally possible, for some P,Q,p,q, one can convert the 
extrinsic homotopy into an intrinsic homotopy of the form 



in which the path parameter t appears linearly. This means that the linear 
algebra to compute C,D G C Nxm and c,d G C N is done just once at the out- 
set, rather than being repeated at each value of t. This can save a significant 
amount of computation and is also simpler to implement. 

This paper is organized as follows. In we review the extrinsic ho- 
motopies formulated in ^3] for intersecting algebraic varieties, and in S J2.1I 
and ^2.21 we show how to convert these to the linear intrinsic form. A com- 
parison of the numerical behavior of the extrinsic homotopies and intrinsic 
homotopies is presented in 



x(u,t) = R (t)u + x p (t) 



(3) 



H(u, t) := H(x(u, t),t) = + x p (t)) = 0, 



(4) 



H(u,t) = f(t(Cu + c) + (1 -t)(Du + d)) = 0, 



(5) 
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1 Extrinsic Diagonal Homotopies 

Let A C C k and B C C k be as in the opening paragraph, having dimensions 
a and b respectively. We have bounds on the dimension of components of 
A PI B as follows. After renaming if necessary, we may assume a > b. The 
largest possible dimension of A fl B is therefore 6, which happens if and 
only if B is contained in A. We can check this possibility using a homotopy 
membership jT2] test to see if a generic point of B is in A. If so, we have 
A C\ B = B and no further computation is needed. Otherwise, we know 
that the largest possible dimension of A fl B is b — 1. On the other hand, 
because the codimension of A fl B is at most the sum of the codimensions of 
the A and B, the smallest possible dimension of any component of A fl B is 
max(a + b — k, 0). For a particular problem, one might have available some 
tighter bounds on dim(yl fl B), and if so, one can take advantage of that 
knowledge in the algorithm to follow. Accordingly, we introduce the symbols 
h* and h as follows: 



max(a + b — k, 0) < ho < min(dim(any component of A fl B)). (7) 
Unless we have other knowledge, we use the defaults h* = b and h = max(a+ 



Instead of working directly in C fc , we find the intersection Af]B by casting 
the problem into (u, v ) G <C h+k and restricting to the diagonal u — v = 0. 
More precisely, the product X := A x B C C fc+fc is an affine variety of 



b>h*> dim(A n B) 



(6) 



b- k,0). 
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dimension a + b, i.e., an irreducible affine algebraic set of dimension a + b. 
The intersection of A and B can be identified, e.g., [21 Ex. 13.15] or [TO! pg. 
122ff], with In A where A is the diagonal of C k+k defined by the system 



S(u, v) 



Ui - v 1 



(8) 



with (u,v) giving the coordinates of C k+k . 

The initial data consists of witness sets for A and B. That is, our data 
for A consists of a generic system La{u) = of a linear equations and the 
deg(A) solutions {a±, . . . , «deg(A)} C C k of the system 



f A {u) 
L A (u) 



0. 



(9) 



and similarly the data for B consists of a generic system Lb{v) = of b linear 
equations and the deg(B) solutions {Pi, . . . , Pdeg(B)} C C m of the system 



fs(v) 
L B (v) 



0. 



(10) 



Remark 1.1 We are not assuming that A and B occur with multiplicity one 
in the solution sets of their respective systems = and Jb(v) = 0. If 

the multiplicity is greater than one, we must use a singular path tracker [15] . 

The extrinsic algorithm can be summarized concisely by introducing a 
bit of matrix notation. First, let 



w 



u 

v 



e C 2fc , (11) 

and introduce a column vector of "slack" variables z G C fc . Also, define the 
k x k projection matrix 



diag (1, . . . , 1,0, 



;i2) 



k-h 
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Left multiplication by Ph picks out the first h rows of its multiplicand and 
right multiplication picks out the first h columns of its multiplier. Note also 
Ph- Similarly, let Pji be the k x k matrix 



that P\ 



diag (0,...,0,1,...,1,0,...,0 




(13) 



which picks out rows (or columns) j + It is useful to note that 

P -4- P = P 

The formulation of the homotopy requires several random matrices as 
follows. First, we choose generic matrices 



(14) 



where #(/U) is the number of functions in the system /a(x) associated to 
component A, and similarly for #(/b). These are used to define 



Mf A (u) 



(15) 



Note that A x B is an irreducible component of the solution set of the system 
^{w) = 0. Next, we choose A a generic (a + b) x k matrix, and let 



A = [ A -A ] e C (a+6)x2fc 
so Aw = A(u — v ) . Finally, we choose generic matrices 

B G C (a+b)xfc , C G C kx2k d G C fcxl . 



(16) 



(17) 



In all these, a matrix with random complex elements will be generic with 
probability one. 

Since the smallest dimensional nonempty component of A H B is of di- 
mension at least max{0, a + b — k}, it follows from [TBI Lemma (3.1)] that 
we can find the irreducible decomposition of A D B by finding the irreducible 
decomposition of Aw = on X = A x B. For this purpose, we consider a 
cascade of homotopies of the form 



£ h {w,z) 



Aw + BP h 2 
z -P h {Cw + d) 



(18) 
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which is well-defined for any integer < h < k. Denoting the entries of 
z as Z\, . . . , Zk, note that the last row of this matrix equation implies that 
(zh+i, ■ ■ ■ , Zk) = 0. The method for generating a witness superset consists of 
solving £h*(w, z) — and then descending sequentially down the cascade to 
solve £j(w, z) = for j = h* — 1, ... , h . 

The rationale behind the cascade is that the linear system P^Cto + d) = 
is a linear slice that cuts out witness points for solution components of 
dimension h. The vector z is a set of slack variables. A solution point 
of £h(w,z) = for which z = is on the slice and thus gives a witness 
point. Solution points with z ^ are not on the slice, and we call these 
"nonsolutions." These become the starting points for the next step of the 
cascade. (We state this more formally below, after giving more details of 
the algorithm.) For each step down the cascade, one more slack variable is 
set to zero and a corresponding hyperplane is removed from the slice. The 
recycling of nonsolutions as starting points for the next step of the cascade 
is valid due to the fact that for j < i, £j(w,z) is just £i(w,z) with certain 
elements of B, C, and d set to zero. This is justified in [To] . 

The following steps of the algorithm still need to be described: 

• how to solve £h*(w, z) = 0, 

• how to descend the cascade, and 

• how to reap the witness points from the solutions at each level of the 
cascade. 



The homotopy to solve £h*(w, z) = is 

T{w) 

A i/i -I- RP 



Aw + BP h *z 
z-P h *(Cw + d) 



L A (u) 
L B (v) 

z 



(19) 



where 7 is a random complex number. At t = 1, solution paths start at the 
deg(A) x deg(i?) nonsingular solutions {(ai, Pi), . . . , {a deg ( A ) , Pde g (B))} C C 2fc 
obtained by combining the witness points for A and B. At t = 0, the solution 
paths terminate at the desired start solutions for £h*(w,z) = 0. In we 
ended the homotopy at £&(u>, z) = 0, but the argument works equally well 
with h* in place of b. 
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The homotopy connecting Si to Sj for j < i is 



Hij(r,w,z) 



T{w) 

Aw + BPiZ 

2-(P i + rP i< )(Cw + d) 



(20) 



where r goes from 1 to along a sufficiently general 1-real-dimensional curve. 
For example, for all but finitely many 7 e C of absolute value 1, r = r + 
7r(l— r) as r goes from 1 to on the real interval suffices. Another possibility, 
relevant in what comes below, is 

r = */(* + 7 (l-*)) (21) 

as * goes from 1 to on the real interval. 

In the cascade of homotopies from [T^] (based on [TT]). we start out with 
the finite set (/, of nonsingular solutions of Si with Zi 7^ 0. Tracking these 
start solutions we end up with a set of solutions Qfj of Sj with Zh — for 
h > j. In |15j . j = i — 1, but the argument there works immediately for any 
j < i. The key points about the set C/f,- is that 

1. the set (/j equals the set of points in Qfj for which Zj 7^ 0; 

2. the set of points Wj C £/f„- for which Zh = for all h < j contains a 
witness point set Wj for the j-dimensional components of the solution 
set of the intersection of A and B. 

We also know that the set of points in Qf. for which z^ = for all h < j 
equals the set of points in yf , for which Zj =0. We wish to set up an intrinsic 
homotopy such that analogs of the above key facts hold true. 



2 Setting Up Intrinsic Homotopies 

The extrinsic homotopies of (|19|) and (|2(Jj) use the variables (w, z) e C 2k x C fe . 
Each has a + b + k linear equations which we wish to eliminate by converting 
to an intrinsic homotopy. The result will be homotopies in intrinsic variables 
y e C 2h ~ a ~ b . Note that 2k — (a + b) is the codimension of A x B in C 2k . 
It is also the sum a + b of the codimension a = k — a of A in C fc and the 
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codimension b = k — b of B in C k . Since this quantity appears frequently in 
the expressions below, we define 

m = 2k-a-b. (22) 

Accordingly, our intrinsic homotopy variables are y G C m . 

2.1 Intrinsic Start Homotopy 

In this section, we replace the extrinsic start homotopy of (JT§j) with one 
having the intrinsic form of (jSJ). Fixing a particular solution 



of 



Wi 



L A (u) 
L B (v) 



Up 



0. 



choose a basis W\ G C 2kxm of the null space N\ of 

= 0. 



L A (u)-L A (0) 
L B {v)-L B (0) 



(23) 



(24) 



(25) 



arising from (J0| and (fTUj) correspond to Ni H 



The solutions (ati, of 
(4xB). 

Fixing a particular solution u> 2 of 

Aw + BP h .(Cu; + d) =0, 
choose a basis W2 G C 2fcxm of the null space N 2 of 

Aw + BP^Cu; = 0. 
We have the intrinsic homotopy with variable y G C m 

T ((1 - r) [tu! + Wiy] + r [ W2 + W 2 y]) = 0. 



(26) 



(27) 



(28) 



Since Aq is transverse to A x B, the (2A; — a — 6)-dimensional affine subspace 
given by 

{ n [ Wl + W lV ] + r 2 [w 2 + W 2 y] I y G C m } (29) 
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is transverse to A x B for all but a finite set of [t 1 ,t 2 ] G P 1 . In particular 
for all but a finite number of 7 G C of absolute value one, with the relation 
between r and t as in (j21j) . the m-dimensional affine subspace given by 

{(1 - r) [ Wl + W x y\ +t[w 2 + W 2 y] | y G C m } (30) 

is transverse to A x B for all t G (0,1]. By genericity in the choices of 
A, B, C, d, this is true for t = also. Thus using the homotopy (|2H|) to track 
the paths starting with the (««, (5j) at t = 1, we get the start solutions of the 
cascade at t = 0. 

In practice it will be convenient to go directly from solutions (ai,/3j) 
of (jUJ) arising from (jHJ) and (fTUj) to or any £j with j < h*. Doing 

this we want to know that the limits of the paths of the intrinsic homotopy 
starting with the solutions (a i; (3j) contain the subset Qj for which Zj ^ 
and a set of points Wj which contains a set of witness points Wj. This is true 
for both the intrinsic and the earlier extrinsic homotopy of ^H]- The reason 
why this is so is that the solutions QjUWj are contained in the set of isolated 
solutions of Sj restricted to Ax B. Therefore by ITol Lemma A.l], there is a 
Zariski open set of t G C such that except for a finite choice of 7 of absolute 
value one in (J2*T|) . Qj U Wj are limits of isolated solutions of the homotopy 
(J2HJ) restricted to Ax B. Since the solutions at t = 1 of the homotopy on 
AxB are the transversal intersection with the m-dimensional affine subspace 
given by Eq.(|30|). it follows that for the t near 1 this is still true. Thus the 
isolated solutions of the homotopy (|28j) for a Zariski open set of the t are 
continuations from solutions (at*, /3j) of arising from @ and (jlOj). and 
in consequence Qj U Wj are contained in limits of isolated solutions of the 
homotopy (j2Hj) restricted to A x B starting at these points. 

The current default is to go directly from solutions (atj, (3j) of ()24j) arising 
from © and (HH to £ h »-i. 

2.2 Intrinsic Cascade Homotopies 

In this section, we convert the extrinsic cascade homotopies of (J20)) into 
intrinsic the form of (JEJ). This must be done a bit more delicately than what 
was done for the start homotopy, because we must preserve the containment 
of Hij inside the parameter space of Si so that we retain the properties stated 
at the end of ^ We do this by deriving an intrinsic homotopy whose path 
is exactly the same as a generic real path from r = 1 to r = in (J20j) . 



10 



We start by eliminating z by substitution from the last block row of (|2(Jj) 
into the middle row. We use the facts that for i > j, PjPj = Pj and 



PjPji = Pji to obtain 



H itj (t,w) :-- 



Aw + B(Pj + rP ii )(Cw + d) 



0. 



(31) 



which, abusing notation, we still call "H%y By similar abuse of notation, we 
use Sh{w) in place of Sh(w, z) after eliminating z from (fTSj). 

Our first observation concerns the existence of a constant particular so- 
lution throughout the cascade. 



Lemma 2.1 The inhomogeneous linear system 



Ik ~Ik 







c 


w = 


-d 



(32) 



has a unique nonzero solution e. 

Proof. The genericity of C implies the invertibility of 



Ik —Ik 
C 



□ 



Notice that this implies that both Ae = and Ce + d = 0, and therefore 
w = e is a solution of 



Aw + B(P j + TPji)(Cw + d) = 

for any i, j, t. 

Let Y/j be the homogeneous linear system 

Y fe :=(A + BP A CW = 0. 



(33) 



(34) 



The following lemma concerning the null space of is crucial for the con- 
version to an intrinsic form. 

Lemma 2.2 For any j and i such that ho < j 1 < % < h* , there exist matrices 
E G C 2kx( - m - i+ ri and F,G G C 2kx ^ such that 

1. [E F] = Null Yi 

2. [E G] = Null Yj 
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3. PjiCF — PjjCG 



h 







where the [i — j) x (i — j) identity matrix Ii-j appears in rows j + 1, . . . , i. 

Proof. We must first establish that Yj and Yj are full row rank a + b so 
that m = 2k — a — b is the correct dimension of their null spaces. Since A 
depends on generic A (see ([Top and B and C are generic, it suffices to show 
that there is at least one choice of A, B, C such that Y^ is full rank for 
ho < h < h*. For a + b < k, it suffices to choose B = 0, C = and choose A 
to make Y^ = [I a +b — I a +b 0]w. For a + b > k, choose A = [J*. 0] T , choose 
B with Ik-a-b in the lower left and C with Ik- a -b in the upper left. Since 
h > k — a — b, this suffices to make Y^ full rank, as one may check by direct 
substitution. 

Next, we establish that Yj and Yj share a null subspace of dimension 
m — i + j. Note that 



Yi = (A + B(P j + P^C)™ = (Y, + BPjiC)w. 



(35) 



The matrix BPjjC is independent of BPjC because the projection matrices 
pick out different rows and columns of generic matrices B and C. Accord- 
ingly, the subspace Null YiRNull Y^ = Null Y^nNull (BP^C). These have 
dimension m and (i — j), respectively, and they meet transversely, so the 
intersection has dimension m — i + j. Let E be any basis for this subspace. 

Now, suppose F completes a basis [E F] for Null Yj. It must be inde- 
pendent of Null (BPjjC), and since B is generic, this implies that PjiCF 
must be full rank. Since P^ zeros out all but rows j + 1, . . . ,i, this implies 
that 



PjiCF 





Q 
o 



(36) 



must have a full-rank (i — j) x (i — j) matrix Q in rows j + 1, . . . , i. Then, 
F = FQ^ 1 completes the basis of Yj while also satisfying Condition 3 of the 
lemma. Similar reasoning shows the existence of G. □ 



Choosing a random 7 G C, we form the linear system 
Wu(t,y) = e+ [ E tF + j{l-t)G]y 



(37) 
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where y G C m . From this, we form the intrinsic homotopy 

H iJ (t,y)=r(W iJ (t,y)) = 0, (38) 

and track y as t goes from 1 to on the real interval. 

The crucial fact behind the equivalence of the intrinsic and extrinsic ho- 
motopies is that the space intrinsically parameterized in ()37)1 is the same for 
appropriate choices of parameters as the space that we extrinsically cut out 
with linear equations before. 

Lemma 2.3 For all but a finite number of 7 G C of absolute value one, it 
follows that for any t G [0, 1] there is a 7^ r G C such that the kernel of the 
linear system 

Aw + B(Pj + rP ji )(Cw + d) = 0. (39) 
on C 2fc is parameterized by Wi j(t,y) where y G C m . 

Proof. This follows immediately for t — and 1 with no restriction on 7 
of absolute value 1 by taking r equal to and 1 respectively. 

Combining this with the dimension of the kernel of (|3T)|) being at least 
m, we conclude that the dimension of the kernel of (J3U|) is exactly m except 
for finitely many ^ 7 G C. In particular, for all but a finite number 7 of 
absolute value 1, the dimension of the kernel of (|39|1 for r and t as in (|2"T|) with 
t G (0, 1) is of dimension m. Since e satisfies both Ae = and Ce + d = 0, 
it is therefore enough to show that for all (t, y) 

(A + B(Pj + rPji)C) [ E tF + j(l-t)G]y = 0. (40) 

Since the columns of E are in Null Yj D Null (BP^C), it is annihilated. 
Since y is arbitrary, we must have 

(A + B(P, + rP^C) [tF + 7(1 - t)G] = 0. (41) 
Since F is in Null Yj and G is in Null Yj, this is the same as 

B ((r - l^PjiCF + r 7 (l - t^^CG) = 0. (42) 
By Condition 3 of Lemma 12.21 this becomes 



((r- l)t + r 7 (l -t))B 





i— _ 
0' 



(43) 
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which equals zero by (|21|1. 



□ 



We rephrase Lemma \2. 31 

Lemma 2.4 For all but a finite number of j G C of absolute value one, it 
follows that for any t G [0, 1], the system 

F(W id {t,y)) = (44) 

on C m is the intrinsic system associated to the system 

Aw + BPiZ 
_ z- (Pj +rP ji ){Cw + L) 

with r = t/(t + 7(1 - t)). 

We define Qi as the set of nonsingular solutions of Hij(l,u) on which 
Pi(Cw + L) is nonzero and which correspond to points of A x B; Qj as the 
set of nonsingular solutions of i?ij(0, uj) on which Pj(Cw + L) is nonzero and 
which correspond to points of A x B; and Qij as the sent of limits obtained 
by tracking Qi from t — 1 to t — using the homotopy Hij(t, uj). 

Theorem 2.5 The subset Wj C Qij on which Pj(Cw + L) is zero contains a 
set of witness points for the j -dimensional components of An B. These wit- 
ness points include deg(Z) distinct points for each irreducible j -dimensional 
component Z of An B. Moreover Qj C Qij. 

Proof. The sets Qi, Qj considered as sets of solutions of the extrinsic systems 
Si, Ej on C 2fc are the same as the sets occurring in the homotopy of jT3] . The 
extrinsic homotopy from jT3] that we discussed in ^is simply a differentiable 
path P parameterized by t G [0, 1] on a complex line I in the parameter space 
of the systems z) joining a general point E; L to a general point Ej of the 
linear subspace of systems of the from Ej(w,z). The only fact about the 
path P used in [12] is that it depends on a choice of 7 G C of absolute 
value 1, which can be chosen, except for a finite number of complex numbers 
of absolute value 1, so that P avoids a certain finite subset B of i. In 
Lemma 12.41 we show that the intrinsic homotopy leads to systems on the 
same complex line £. What changed is that the path P' on I is not linearly 



(45) 
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related to the original path P. But since the path P' depends on a choice 
of 7 G C of absolute value 1, which can still be chosen, except for a finite 
number of complex numbers of absolute value 1, so that P' avoids the finite 
subset B of £, the same conclusions of still hold. □ 

2.3 Algorithm Summary 

The homotopy algorithm to intersect two positive dimensional varieties in 
intrinsic coordinates is described below. After the initialization, there are 
three stages. First is the homotopy to start the cascade, followed by the 
homotopy to find a witness sets for the top dimensional part of A D B. 
Thirdly, all lower dimensional parts of A D B are computed in a loop from 
6 — 2 down to h . The second and third stage are separate because we 
can avoid a coordinate transformation. Also, in many cases - such as the 
important application of the intersection with a hypersurface - the loop will 
never be executed. 

Some subroutines used in the algorithm below are just implementations of 
one formula in the paper, e.g.: Combine implements (II 5|) . Next we describe 
briefly the other subroutines. 

The linear algebra operations to deal with solutions in intrinsic coordi- 
nates are provided in the subroutines Start_Plane, Project, Initialize, 
Basis, and Transform. Given the equations for La and Lb, Start_Plane 
first computes a basis for the null space of L A ^(Q) and Lg l (0) before doubling 
the coordinates into a corresponding basis in C 2fc . After orthonormalization 
of the basis, Project computes the intrinsic coordinates for the product of 
the given witness sets of A and B. The subroutine Initialize first generates 
the random matrices A, B, C, and d before computing the e of Lemma \2. II 
In addition, Initialize returns the operator Y, which returns for any h the 
corresponding Y/j of ([34)1 . Lemma l2~2l is implemented by Basis, while Trans- 
form converts the coordinates for the solutions from one basis into another. 

The path tracking is done by the procedure Track. On input are the 
homotopy and start solutions. Except from the set up of the homotopy in 
intrinsic coordinates, one can implement Track along the lines of general 
path following methods, see [I], [HI El, or jHj. 

The subroutine Filter takes on input the witness sets W for higher di- 
mensional components and the list Z. On return is W, augmented with a 
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witness set for the solution set at the current dimension, and a filtered list 
Z of nonsolutions. The list Z given to Filter may contain points on higher 
dimensional solution sets. To remove such points, a homotopy membership 
test as proposed in [TJ] can be applied. Recently, an interesting alternative 
was proposed by Li and Zeng in Bj. The nonsolutions serve as start solutions 
in the cascade to find witness sets for the lower dimensional solution sets. If 
Z becomes empty after Filter, the algorithm terminates. 



Algorithm 2.6 Intersecting two ] 
Input: k, a, b (a > b); 

fA(u) = 0J B (v) = 0; 
L A (u) = 0,L B (v) = 0; c 

w A ,w B . 

Output: J-{x) = 0; 

L = [Lh , • • • , £&-i]; 
W=[W h0 ,...,W b - 1 ]. 

T := Combine(/ A , f B ); 
S := Start_Plane(L A ,L B ); 
Z :=Project(W A x W B ,S); 
[Y, e] := Initialize^, a, b); 
[E, F, G] :=Basis(Y 6 ,Y 6 _ 1 ); 
W(t,y) := [tS+(l-t)[e+[E F]]; 

Z:=Tr a ck(F(W(t,y)),Z); 
Z := Track^, [E,F,G],Z); 
[W 6 _!,Z] :=Filter(W,Z); 
ho := max(a + b — k, 0); 
for j from 6 — 2 down to ho do 

[E, F, G) :=Basis(Y J+1 ,Y J ); 

Z := Transform^, [E, F]); 

Z := Track^, [E,F,G],Z); 

[Wj,Z] := Filter (W,Z); 
end for. 



'ositive Dimensional Varieties A and B. 
dim(A) = a, dim(S) = b, A, B C C k 
polynomial systems in u,v G C k 
im(L/(0)) = k-a, dirn^^O)) = k-b 
solutions in witness sets for A and B 
system combined from f A , f B in x £ C k 
list of linear spaces, dim(L~ 1 (0)) = i 
solutions Wi in i-dim witness sets 

combine systems f A and f B as in (|15j) 
basis for plane defining Wa x VVb 
solutions to start the cascade 
linear space Aw + BP/ l C(w + d) = 
basis for Null Y b and Null Y b _i 
deform start plane S into [E F] 
with t using formula (}2"Tj) 
homotopy to start the cascade 
find top dimensional component 
keep witness sets and nonsolutions 
minimal dim(A D B) 
compute witness set at dimension j 
W(t,y) = e+[E tF + 7 (l-t)G]y 
coordinates into new basis [E F] 
homotopy FiWj+ijit^y)) = 
keep witness sets and nonsolutions 
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3 Numerical Experiments 



The algorithms in this paper have been implemented and tested with PHC- 
pack PHI- To compare with our implementation in extrinsic coordinates, we 
use the same examples as in jT^j. All computations were done on a 2.4 Ghz 
Linux machine. 

(1) An Example from Calculus. In this example, we intersect a cylin- 

der A with a sphere B. More precisely, A = { (x,y,z) \ x 2 + y 2 — 1 = } 
and B = {(x,y,z) \ (x + 0.5) 2 + y 2 + z 2 — 1 = }. The intersection 
A fl B is a curve of degree four. Since k = 3, a = 2, and 6 = 2: h = 1, 
so there are only two homotopies, each defining four solution paths. 

(2) An Illustration of the Cascade. In this example we need to execute 

the cascade to find the point of intersection. We consider the compo- 
nents A = { x = 0,y = } and B — { z — 0,w — } &s solution sets of 
the same system f(x,y,z,w) = [xz,xw,yz,yw] T = 0. We have k = 4, 
a = 2, and 6 = 2. 

(3) Adding an Extra Leg to a Moving Platform. In this example we 

cut a hypersurface A in C 8 with a curve B, i.e.: a = 7 and 6 = 1. The 
application concerns a Griffis-Dufty platform [S] (analyzed by Husty 
and Karger in [5] and subsequently in jHj) where A fl B can be inter- 
preted as adding a seventh leg to the platform so it no longer moves. 
As deg(A) = 2 and deg(5) = 28 (ignoring the mechanically irrelevant 
components), there are 56 paths to trace, by two homotopies. 

In the Table below we list all important dimensions of the three example 
applications. A summary of the execution times is reported in Table 121 

In these numerical experiments, we save about half of the computational 
time when working in intrinsic coordinates. Comparing the number of vari- 
ables of the original extrinsic method, M = 3k — a for the examples tested, 
with the number for the intrinsic method, m = 2k — deg(A) — deg(B), we 
have in these experiments 3k — a = 7, 10, 17 variables reduced to 2,4,8, or 
more than half. Since the cost of linear solving is 0(n 3 ), this implies about 
a eight-fold reduction in the cost of that portion of the code. Linear solving 
can be a significant portion of the total cost, as it is used in Newton's method 
for tracking the homotopy paths. The experimental results suggest that this 
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example 


dimensions and degrees of A and B 






deg(A) 




k 


dim(A) 


deg(A) 


dim(S) 


deg(5) 


m 


M 


x deg(.B) 


(1) 


3 


2 


2 


2 


2 


2 


7 


4 


(2) 


4 


2 


1 


2 


1 


4 


10 


1 


(3) 


8 


7 


2 


1 


28 


8 


17 


56 



Table 1: Dimension and degrees of the two irreducible sets A and B for 
the three examples, followed by # variables m = 2k — dim(A) — dim(S), 
M = 3k — a (which is the ^variables in the extrinsic homotopy), and number 
of paths deg(A) x deg(B) at the start of the cascade. 





Homotopies 


Total CPU Time 





1 


2 


intrinsic 


extrinsic 


(1) 


0.03 


0.01 




0.04 


0.07 


(2) 


0.01 


0.02 


0.01 


0.04 


0.11 


(3) 


9.90 


5.94 




15.84 


34.70 



Table 2: Timings in CPU user seconds on 2.4Ghz Linux machine. The 
second column concerns the homotopy to start the cascade, in the third 
column are the timings for the top dimensional components, followed by the 
eventual next homotopy in the cascade. 



was accounting for about half of the total cost in the extrinsic method, but 
accounts for a much less significant fraction of the computational cost of 
the intrinsic method. The other 50% or so of the cost remains, which is 
attributable to function evaluation, data transfer, and other overhead. The 
cost of function evaluation can vary dramatically from one polynomial sys- 
tem to another, so we cannot definitively expect the same percentage savings 
for all systems, but we can say that the intrinsic formulation seems to give a 
substantial reduction in computational time. 
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